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ABSTRACT 

The high multiplicity of massive stars in dense, young clusters is established early in their evolu- 
tion. The mechanism behind this remains unresolved. Recent results suggest that massive protostars 
may capture companions through disk interactions with much higher efficiency than their solar mass 
counterparts. However, this conclusion is based on analytic determinations of capture rates and esti- 
mates of the robustness of the resulting binaries. We present the results of coupled n-body and SPH 
simulations of star-disk encounters to further test the idea that disk-captured binaries contribute to 
the observed multiplicity of massive stars. 

Subject headings: binaries: general — circumstellar matter — stars: formation 



1. INTRODUCTION 

Massive stars such as those of the Trapezium in 
the Orion Nebula cluster (ONC) have a hig h mul- 
tiplic ity compared to l ow m as s stars (iMason et all 
1991: iPreibisch et all fl999l: iStahler et all 120001 : 
Garcia fc Mermilliodl 1200 The formation of a 

multiple system in the early life of a cluster can occur 
through fragmentation of the prestellar material, or by 
capture. Capture can occur via multi-body interactions 
with other cluster members, or through disk interactions 
in the protostellar phase. There is growing observational 
evidence for massi ve, embedded disks surrounding 
massive protostars (|Cesaroni et all 120071 . for a recent 
review). Simulations of massive star formation suggest 
that the masses of the disks may build up to values ~ 
30% of the central star's mass b efore global instabilitie s 
trigger a sudden accretion event (|Krumholz et all 20071 
Fragmentat ion of such massive disk s into companions 
is possible ()Kratter fe Matznerl ^0061 but in this work 
we consider capture of a lower mass cluster member by 
a massive star-disk system , contin uing the analysis pre- 
sented in IMoeckel fe Ballvl (|2006l) and IMoeckel fe Ballvl 
(127)071 hereafter MB07). 

Disk assisted capture in low mass systems has been , 
studied using analytic methods (|Clarke fc Pringldll991h 
and smoothed particle hydro dynamics (SPH) codes 
(|Hellerlll995t iBoffin et allll998l l The capture rates from 
these studies are too low to be a significant contributer to 
the binarity of roughly solar mass stars; for 20 M Q stars 
with a 2 Mq disk in a Trapezium- like environment, MB07 
show that the capture rate is high enough to account 
for 50% binarity after 1 Myr. MB07 estimated analyti- 
cally the likelihood of survival for these binaries, and con- 
cluded that higher mass captured companions are more 
likely to survive, while lower mass companions would be 
preferentially ionized by encounters with other cluster 
members. However, in calculating capture rates and es- 
timating survival fractions, many averages and integrals 
are taken. It is desirable to further test these results in 
an unaveraged, more 'realistic' setting. 
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In this paper we describe the results of n-body simula- 
tions of a cluster similar to the ONC, in which we include 
the dissipative effects of passages through a circumstellar 
disk surrounding a central 20 M Q star. This test demon- 
strates that the rates derived by MB07 are reasonable in 
a cluster setting, and that the most likely outcome of en- 
counters between the captured binary and other cluster 
members leaves the central star in a binary. 

2. SIMULATIONS 

In this study we combine the SPH results of MB07 
with n-body simulations of a cluster sim ilar to the ONC. 
This i s similar in spirit to the work of IScallv fc Clarkd 
(120011) and Pfalzner et al. dPfalzner et al.l l2006l: iPfalznerl 
200a |Pfarzner fc Olczakll2007D . in which n-body simula- 
tions of a cluster are used to determine encounter fre- 
quency and parameters, which are then used to study 
the effect of encounters on disks. This work is differ- 
ent in that the results of the close encounters are in- 
cluded in the simulation as it is running, similar to 
iMcDonald fc Clarkel ([1995), who used an analytic pre- 
scription for disk-mediated binary formation in simula- 
tions of a cluster with 10 stars. 

The sim ulations are p erformed using Aarseth's code 
NBODY6 (lAarsethl I2000H The cluster is set up as a 
King model (e.g. iBinnev fc Tremaind Il987t ) with Wo — 
9, core radius r c = 0.2 pc, central density n — 2.0 x 10 4 
pc -3 , and central velocity dispersion <jq = 2.19 km 
These parameters are similar to those of the ONC 



(Hillenbra nd fc Hartmann| [l998). and with a cut-off ra- 
dius of 2.0 pc yield 4725 sta r s in t he simulation. The 
IMF used is that of iKroupal (|2002ft in the range 0.3 - 
9.0 M© . Each star is randomly placed according to the 
King model density distribution, with a random veloc- 
ity appropriate to its radius. In addition, we place a 20 
Mq star at the center of the simulation with a random 
velocity. 1000 sets of initial conditions were generated; 
each is run with and without the effects of disk dissipa- 
tion. In the dissipation runs, the only star with a disk 
is the central star. This is an artificial situation; while 
the most massive disk will dominate an interaction, the 
presence of disks around all stars in the cluster would 
increase the effect of disk encounters. 
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Fig. 1. — Energy change as a f unction of remnant disk mass 
for three of the binaries simulated in Mocckel & Bally (2006). The 
outlying point for i = 180° shows that the true effect of disk dis- 
ruption on repeated encounters is more complicated than a simple 
mass scaling. However, for most encounters, a scaling between the 
two bold lines is reasonable. 

MB07 used a modified version of the publicly released 
SPH code GADGET-2 (|Srjringell[200^ to simulate en- 
counters between the massive star-disk system and im- 
pactors with masses in the range 0.3 - 9.0 M Q , perias- 
tra in the range 50 - 550 AU, and inclinations from 0° 
- 180°. The disk radius is 500 AU. Interpolation of the 
data from these simulations gives, for any impactor mass, 
periastron, and inclination angle, a change in orbital en- 
ergy and change in disk mass associated with the en- 
counter. We have modified NBODY6 to detect encoun- 
ters between the central star and other cluster members. 
At the encounter periastron, we determine the orbital pa- 
rameters, and then modify the velocity of both encounter 
partners in a momentum conserving fashion so that the 
change in orbital energy is equal to that found via in- 
terpolation of the data from MB07. The change in disk 
mass due to the encounter is also tracked. Because we 
only have data for periastra > 50 AU, encounters with 
smaller periastra are assumed to occur at 50 AU. 

When the orbit of a relatively massive impactor is 
coplanar to the disk, MB07 find that accretion and disk 
capture can increase the impactor's mass by up to 10%. 
However, the change in orbital energy during these pas- 
sages is an order of magnitude greater than can be ac- 
counted for by accretion drag; the dominant contribution 
to the energy change is the change in velocity due to disk 
interactions. We make the simplification that the change 
in orbital energy is due entirely to a change in the rel- 
ative velocity of the two stars. Under this assumption, 
at periastron the velocity kick <5v on the impactor for a 
given orbital energy change AE is given by 

Jv = l\-2v+ Jav 2 + ,, 8A£ .„. y, (1) 
2 ^ y m(l + m/M) J ' y ' 

where v is the pre-kick velocity, m is the impactor mass, 
and M is the primary mass. The change in velocity for 
the primary is 5V = —(m/M)5v. The change in the 
total energy due to these kicks is tracked and included 
in energy checks. 



Fig. 2. — The distribution of encounter frequency for runs with no 
disk dissipation. Also plotted are the best fit Poisson distribution 
to the realizations in which the central star remains in the cluster 
core (solid), and the expected value from equation \3\ (dotted). The 
error bars are single sided 1-er confidence limits (Gchrcls 1986) 

In order to account for the change in disk mass, we 
scale the change in energy for an encounter according to 

AE(m d ,i,mi,r p ) = AE Q (i,mi,r p ) ( — ^- ) . (2) 

\m d0 J 

Here AE is the orbital energy change from an encounter 
with disk mass nid , impactor mass rrii , inclination angle i , 
and periastron r p . AEq is the energy change for the same 
parameters with the disk at its original mass m^o, and is 
found from interpolation of the data in MB07. The index 
n we take to be 1 or 2, and scales the proportionality 
of the the energy change with the remaining disk mass 
fraction. 

The scaling of energy ch ange with d isk mass is not 
completely straightforward. Heller (1995), simulating en- 
counters with disks of different masses, found that the 
change in energy scales r oughly linearly wi t h the disk 
mass, in which case n = 1. IMoeckel fc Bailvl ((2006) sim- 
ulated repeated encounters between a captured impactor 
and the same remnant disk. Analysis of that data (Fig. 
[T]) shows that the scaling of energy change and disk mass 
is more like n = 2. However, this also takes into ac- 
count the changing radius of the disk, an effect which 
we do not include in this study. Since the proper scaling 
depends on details of the encounters that would require 
full SPH simulations of each close passage, we instead 
run all simulations with both scalings and compare the 
results. 

Our scheme prese rves the binary pe r iastro n separation, 
which is shown by IMoeckel fc Ballyj (|2006[ ) to decrease 
during repeated retrograde encounters, and increase with 
prograde encounters. In the most extreme (and rarest) 
case, an in-plane retrograde passage, the periastron sepa- 
ration decreases by ~ 25% after the first encounter. The 
data in MB07 show that the change in energy scales com- 
parably to the periastron radius; thus we would expect 
an error < 25% due to our artifically fixed periastron. 
This is similar in magnitude to the uncertainty in our 
mass scaling, which as shown below has a negligible im- 
pact on our results. 
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Each case is run for 0.5 Myr; we limit the analysis 
to this time for three reasons. After this time the disks 
are mostly destroyed. In our simulations this destruction 
is by encounters, while in reality the additional effect of 
photo-evaporation will contribute. Mass segregation also 
begins to take effect after this time, which is an added 
complication in the analysis of the results. Finally, the 
multiplicity of s tars is established early in their evolution 
(Mathicu 1994), and a mechanism for binary formation 
should work on short enough timescales to reflect that. 

3. RESULTS 
3.1. Encounter Rates 

The calculation of binary capture rates (Clarke & 
Pringle 1991; Heller 1995; Boffin et al. 1998; MB07) is 
largely similar to the estimation of collision or encounter 
rates. We begin by comparing the standard encounter 
timescale calculation to the simulations. 

The encounter rate for a star of mass M in a cluster 
with number density n, velocity dispersion c, and en- 
counter radius r is given by 



7 = 4v / 7iV 



1 



GM 



(3) 



Here M. is M + fh, with fh the average stellar mass in 
the cluster (Binney & Tremaine 1987; MB07 for general 
masses). For the parameters of our simulated clusters, 
the encounter rate is 7 = 1.02 x 10~ 5 yr -1 , with 5.1 
encounters expected over 0.5 Myr. 

Plotted in figure [2] is the distribution of the number 
of unique cluster members encountered by the central 
star, without disk dissipation. Multiple encounters with 
the same star, for instance in a binary, are counted only 
once. Because some of the central stars have a high ini- 
tial velocity and escape to the cluster outskirts, there is 
an excess of low encounter-number runs. Therefore the 
distribution for all runs is plotted, as well as only those in 
which the central star remained within the cluster core. 
One would expect that the distribution for stars in the 
same environment would be Poisson; the data shows oth- 
erwise. Because the number density and velocity disper- 
sion are radially dependent, central stars that move to 
larger radii are exposed to a much different environment 
than those that remain near the cluster center. By limit- 
ing our analysis to the cluster core, where the properties 
are closest to those used in equation |3l the best fit (with 
mean value 5.78) is in reasonable agreement with the 
theoretical value of 5.1. 

Equation[3]is averaged over the mass function. Because 
of the dependence of the encounter rate on the mass of 
the stars, encounters with more massive cluster members 
are more frequent. The mass function of encounter part- 
ners is well fit by a single power-law mass function of the 
form £(m) cx mT a with a = 0.92, while a Salpeter mass 
function has a = 2.35. 

3.2. Binary Fraction and Mass Function 

Of greatest interest is the fraction of the massive stars 
at a given time that are in a multiple system. MB07 
calculate that for a cluster with the parameters of our 
models, the binary formation rate is T = 0.6 Myr -1 
(Figure 3 in MB07). We compare this to our simula- 
tions as follows. For each simulation, a list of stars that 



O n = 1 
A a = S 
□ no disks 
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Fig. 3. — The number of runs with the central star in a binary, 
Nsit), as a fraction of the total number of runs Nt for the three 
simulation series. Cases with both n = 1 and n = 2 in equation [2] 
are at ~ 30% after 0.5 Myr. 



encounter the central star is generated, with the times of 
their encounters. If the same star is involved in consec- 
utive encounters, a binary has formed, and we track the 
following events. 

Ionization: The next two encounters are with differ- 
ent stars. Even if the intruder and the original binary 
partner form a binary, the central star is no longer in a 
multiple system. Exchange: The next two encounters are 
with the same star, which is different from the initial bi- 
nary. Flyby: The next encounter is with a different star, 
but the one after that is with the initial binary partner. 
In the latter two cases, the central star is considered to be 
in a binary throughout. In this scheme the possibility ex- 
ists that random, consecutive encounters with the same 
star could contaminate the binary statistics. In practice, 
the number of binaries with semi-major axes larger than 
the average inter-stellar distance in the cluster core is on 
the order of 1%, and we consider the detected binaries 
to be true binaries. 

Plotted in figure [3] are the number of central stars in 
binaries as a function of time, for each of the three series 
of simulations. Considering first the series with no disks, 
we see an initial rise in the binary fraction, followed by 
a leveling off to ~ 12% at 0.5 Myr, as binary creation 
through random dynamical processes is balanced by ion- 
ization. For the two series with disk dissipation, there is a 
steady rise up to a value of approximately 30%, in good 
agreement with the calculations in MB07. The binary 
fraction does not appear to depend on the specifics of 
the energy-change disk-mass scaling; the early passages, 
when the disk still has nearly its original mass, account 
for the increased binarity. 

Since massive impactors are preferentially captured by 
the disk and more likely to remain bound during ex- 
changes and flybys, the mass function of binary partners 
at 0.5 Myr is flatter than that of the encounter partners. 
With no disks, the best fit single power-law mass func- 
tion has a — 0.67. The disk case with n = 1 has a 
0. If), and for n — 2 we have a = 0.12. 
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TABLE 1 
Binary fate statistics at 0.5 Myr 



Series 


N B 


Ne 


Ni 


Np 


N S (N E - 


f~ Nf)/Ni 


no disk 


387 


44 


224 


431 


119 


2.12 


n = 1 


708 


186 


217 


1619 


305 


8.32 


n = 2 


698 


169 


238 


1482 


291 


6.94 



3.3. Binary Robustness 

MB07 estimated the survival probabilities of capture- 
formed binaries, finding that massive companions are 
more likely to survive encounters with other cluster mem- 
bers, and that lower mass companions are likely to be ion- 
ized. The binary fractions found here are in agreement 
with calculations that don't include ionization effects; in 
order to explain this we turn to the question of ionization 
versus exchange and flybys. 

Shown in table [1] are Nb, the number of binaries 
formed, Ne, the number of exchanges, Ni, the number 
of ionizations, and Np, the number of flybys that occur 
in each series of 1000 simulations. Nb includes all unique 
binary pairings, so that an exchange contributes twice. 
Thus the number of surviving binaries Ns is given by 
Nb — Ne — Nj. The ratio of exchanges and flybys to 
ionizations, (Ne + Ne)/Nj, is indicative of the survival 
chances of a binary in each series. For the diskless case, 
n = 1, and n = 2 this ratio is 2.12, 8.32, and 6.94 respec- 
tively. Encounters between binaries and other cluster 
members are much more likely to end with the central 
star in a binary for the simulations with dissipation com- 
pared to the diskless case, and suggest that the binaries 
formed via this capture mechanism are more robust than 
indicated by the simple estimations of MB07. The in- 
creased disk dissipation in the n = 1 series yields slightly 
harder binaries, but the total binary fraction is not sig- 
nificantly affected by the scaling. 

4. DISCUSSION 

The simulations presented here are intended to test and 
verify the capture rates calculated in MB07. The conclu- 
sions of that work are largely upheld; encounters occur 
at approximately the expected frequency, and binaries 
are captured at a rate consistent with the analytical es- 
timates. In addition, once a binary is formed it is less 



likely to be destroyed by ionization than the estimates 
in MB07 suggest. Encounters between the captured bi- 
naries and intruding cluster stars are far more likely to 
result in a flyby or exchange than in an ionization, leav- 
ing the central star in a binary. 

The capture rates for the central star-disk system are 
not high enough to fully account for the high multiplicity 
observed in massive, cluster-bound stars. Our simula- 
tions produce ~ 30% binarity at 0.5 Myr, the time when 
the disks are mostly destroyed by encounters or photo- 
evaporation, an effect not modeled here. Recent work 
(jKrumholz et al.ll2007h shows that during the formation 
of a massive star, material moves through the protostel- 
lar disk in sporadic, massive accretion events, between 
which the disk builds up to large masses. A disk that is 
~ 30 — 50% of the central star's mass, instead of the 10% 
used here, could increase the capture rates by a factor of 
several. Additionally, continued accretion onto the sys- 
tem or the presence of disks around all the stars could 
increase the capture rates. It is worth noting that our 
simulations here are tailored to ONC-like systems. Since 
the capture rate is linear with stellar density and drops 
with higher velocity dispersion, changing the cluster pa- 
rameters will affect the results. 

The effect of encounters on accretion processes in mas- 
sive star formation is unclear. It is possible that the de- 
struction of the disk by repeated passages could truncate 
accretion at the time of binary formation. Alternatively, 
accretion could resume onto the binary and tighten the 
orbit, leading to a massive, close binary system. The 
frequency of the encounters modeled here is high enough 
that such a situation warrants further investigation. 

As concluded in MB07, this rate can not be ignored. 
However, additional binary formation mechanisms must 
be employed to explain the observed multiplicity of mas- 
sive stars. In the Trapez ium the massive stars h ave, on 
average, 1.5 companions (jZinnecker fc Batdl2002t ). bmce 
disks are effectively destroyed during binary capture, this 
is a process that can only account for a single compan- 
ion, unless further accretion creates a new circumstellar 
or circumbinary disk. 
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